Rhythmogenic neuronal networks, pacemakers, and k-cores 
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Neuronal networks are controlled by a combination of the dynamics of individual neurons and the 
connectivity of the network that links them together. We study a minimal model of the preBotzinger 
complex, a small neuronal network that controls the breathing rhythm of mammals through periodic 
firing bursts. We show that the properties of a such a randomly connected network of identical exci- 
tatory neurons are fundamentally different from those of uniformly connected neuronal networks as 
described by mean-field theory. We show that (i) the connectivity properties of the networks deter- 
mines the location of emergent pacemakers that trigger the firing bursts and (ii) that the collective 
desensitization that terminates the firing bursts is determined again by the network connectivity, 
through k-core clusters of neurons. 



PACS numbers: 87.19.L-, 87.10.Ij, 05.45.Xt 

A neuronal network is a group of interconnected neu- 
rons functioning as a circuit. Each neuron receives elec- 
trical signals from a collection of tree-like dendrites, con- 
nected via synaptic junctions to the branched output ter- 
minals of other neurons. The neuron responds, based on 
some function of its input signals, by either doing noth- 
ing or by "firing," i.e., by producing an action-potential 
output pulse that is received by other neurons [T]. In a 
network of excitatory "integrating" neurons, the electri- 
cal potential of the cell body of a neuron always increases 
by an amount when the cell receives a voltage input 
pulse. The potential of the cell effectively integrates the 
signals from other neuronal outputs. The firing proba- 
bility of a neuron depends sensitively on the electrical 
potential of the cell, leading to threshold behavior in 
which the neuron can be considered to be either in a 
quiescent state characterized by sporadic firing if the cell 
potential is large and negative ( "hyperpolarized" ) , or an 
activated state at higher potentials ( "depolarized" ) , with 
more than an order of magnitude increase in firing rate 
over the quiescent state. 

A classical example of an integrating neuronal network 
is the preBotzinger Complex (pBC) of about 10^ neurons 
located in the brain stem [21 [3]. In this network, which 
collectively produces a rhythmic voltage signal that sets 
the timing of inspiration in mammals under resting con- 
ditions, each neuron is connected on average to one-sixth 
of the other neurons. The period of the current bursts 
is on the order of a second, which is about 10^ times 
longer than the time scale associated with repeated firing 
by activated individual neurons. The slow modulation 
is believed to be due to calcium-mediated "adaptation." 
With each input pulse, the dentritic calcium concentra- 
tion increases by an amount AC. The increase in cal- 
cium concentration leads to an increase in the leakage 
conductance between the dendrites and the surrounding 
medium, making the neuron's somatic potential insen- 



sitive to incoming action potentials. When the somatic 
neuron potential drops below the threshold, it stops fir- 
ing. After a recovery period during which the dendritic 
calcium relaxes back to its equilibrium value, the neuron 
once again begins integrating input signals. 

Two different mechanisms have been proposed to ex- 
plain the synchronization of the firing of the different 
neurons of the pBC. Neurons that in isolation can oscil- 
late autonomously between cycles of firing and quiescence 
are known as pacemakers @]. In the individual pace- 
maker hypothesis, it is assumed that the pBC neurons 
are slaved to a small number of pacemakers believed to be 
present in the pBC. In the emergent pacemaker hypothe- 
sis (EPH) , it is assumed that the oscillation is a collective 
property of a large group of neurons that would not oscil- 
late in isolation [S]. True pacemaker neurons, if present, 
only provide a back-up function and are not essential for 
the oscillation. The oscillations are in this case expected 
to disappear if the number of neurons drops below some 
threshold value. Indeed, when more than 80% of these 
pBC neurons are destroyed in an in vitro experiment, the 
firing sequence changes from periodic oscillation to an in- 
creasingly complex pattern. This also happens when the 
excitability of the neurons is increased [Ej. 

In this letter we show that the ability of a non-uniform 
neuronal network to collectively generate an oscillation, 
as assumed in the EPH, depends on general properties 
of the network connectivity, independent of the details of 
the model of neuron dynamics; moreover, these proper- 
ties can be analyzed in terms of simple graph-theoretic 
methods [7, '5]. Specifically, we show that: (i) a network 
of identical but randomly connected neurons supports 
periodic synchronized bursts triggered by those neurons 
that are linked to a maximum number of other network 
neurons through a minimum number of links, and (ii) 
for highly excitable networks, the minimum number of 
neurons required for rhythmogenesis is determined a se- 
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quence of Magical Numbers, which are a function solely 
of the adjacency matrix. 

We demonstrate these claims using a simple model for 
an excitatory neuronal network regulated by calcium- 
based adaptation. Each neuron is represented by two 
dynamical variables: the somatic potential Vi(t) and the 
calcium concentration Ci{t) of the i**^ neuron. The N 
neurons fire according to the 2N coupled non-linear rate 
equations: 



dt Tv 

da _ 1 
dt TC 



(Kq - V,) + AV{Q) J2 M..,P{V,) (1) 



(Ceq-a) + AC^M,,P(y,) (2) 



Vcq and Ccq are, respectively, the resting potential and 
the equilibrium calcium concentration of a neuron with 
Tv (10ms) and tc (500ms) Hj the respective equilibration 
times (the calcium concentration is thus the slow vari- 
able) . Calcium-mediated adaptation is allowed for by as- 
suming that AT^(C) drops rapidly when the calcium con- 
centration C exceeds a threshold C* . The time-sequence 
of firing events generated by a neuron is assumed to be 
a Poisson process with a voltage-dependent mean firing 
rate P{V). For P{V) we will assume that if V exceeds 
the threshold V* then P{V) increases from a basal rate 
of about five spikes per second to a high rate of about 
seventy-five spikes per second. Finally, the entries of the 
adjacency matrix Mij are equal to one if the output of 
J*'' neuron is an input to neuron i, and zero otherwise. 

We start with the very simplest case of a homogeneous 
network where every neuron is linked to every other neu- 
ron in both directions: Mij = 1 for all z,j (known as a 
"clique" ) . If the initial potentials and calcium concen- 
tration also are the same for all neurons, then the 2N 
rate equations reduce to a single pair that describes all 
neurons: 
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(T4q -V) + NAV{C)P{V) (3) 
(Ceq-C) + 7VACF(F), (4) 



which can be analyzed by the standard methods of dy- 
namical systems [TD] . The pair of equations can also be 
viewed as a "mean-field" approximation for more com- 
plex networks [II] . The resulting dynamical phase dia- 
gram is shown in the upper panel of Fig. [l] 

In the parameter regime marked SO ("stable oscilla- 
tion" ) , the potential and calcium concentrations of the 
neurons undergo a stable limit-cycle oscillation. For 
lower values of the input voltage jump at zero calcium 
concentration AT^(C = 0), corresponding to weakly ex- 
citable neurons, the period of the oscillation increases 
and then diverges as the number N of neurons is reduced 
due to the appearance of an infinite-period saddle-node 
bifurcation [T^]. A line of these bifurcations separates 
the SO phase from a quiescent phase, marked Q, where 
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FIG. 1: (color online) (A) Phase diagram of a homogeneous 
TV-neuron network with every neuron linked to every other 
neuron. The horizontal axis is the maximum voltage jump of 
a cell following an action-potential input pulse. The red line is 
the stability limit of a low-activity fixed point of Eqs. ([3|,(|4| 
(Q phase) while the blue line is the stability limit of a high- 
activity fixed point (HA phase). Cooperative limit-cycle os- 
cillations are fully stable only in the region above the blue 
and red lines (SO phase). (B) Phase diagram of an inhomo- 
geneous, random A'^-neuron network with, on average, each 
neuron linked to A/6 other neurons. In the section labeled 
HA deterministic chaos, period-doubling and intermittency 
is encountered. The dashed lines mark a sequence of Magi- 
cal Numbers A^ determined by the adjacency matrix of the 
network. 



all neurons are permanently in a state of low activity. 
For higher values of AF(C = 0), corresponding to highly 
excitable neurons, the unstable fixed point at the cen- 
ter of the limit cycle becomes stable as N is reduced. 
In the part of the phase diagram where this happens, 
marked HA, the neurons are permanently in a state of 
high activity. This mean-field HA phase does not show 
the complex firing pattern reported experimentally when 
the excitability was increased [6]. 

In actuality, each neuron of the pBC is believed to be 
linked to (1/6)"^ of the other neurons so the pBC network 
is not a clique. To describe this, we use an Erdos-Renyi 
random adjacency matrix |I31 II4j . assigning zeros and 
ones as the entries of My with probabilities 5/6 and 1/6, 
respectively. Solution of the coupled rate equations on a 
single random graph produces the phase-diagram shown 
in the lower panel of Fig. [l] The heterogeneity of the 
network does not destroy its ability to produce robust, 
synchronized stable oscillations, though note that the SO 
section of the phase diagram has been reduced in area as 
compared to the mean-field case. 

Unlike the mean-field case, the firing pattern now 



3 




1.8 2.0 2.2 2.4 2.6 

time (s) 



^ 
> 




10 15 20 25 
time (s) 

FIG. 2: (color online) (A) Periodic potential oscillations (SO 
phase) for a random connectivity matrix. Different colors 
refer to the different neurons. In the low-activity part of the 
cycle, the potential of all neurons is below the firing threshold 
(— 55mV). The potential of a limited number of neurons (e.g. 
blue) rises significantly more quickly at the initiation of a 
burst. When one of these "pacemaker" neurons crosses the 
threshold, it triggers a voltage avalanche that spreads over the 
whole network. (B) Time-dependent potentials in the high- 
activity phase. Multiple collective potential bursts alternate 
with an incoherent, chaotic state. Note the different time 
scales in panels A and B. 



varies greatly from one neuron to the next. Superim- 
posing the firing patterns of different neurons reveals an 
important feature (see Fig. |2j\). In the low-activity part 
of the cycle, the potentials of all neurons are below V* , 
but they rise more rapidly for a sub-population; these 
reach the firing threshold first. Their increased firing rate 
pushes sub-threshold neurons linked to them past the fir- 
ing threshold as well. A chain reaction spreads over the 
network until all neurons are above threshold. Note that 
the least excitable neurons that crossed the firing thresh- 
old latest remain active over a longer period of time, pro- 
ducing a highly asymmetric pulse shape. Even though in 
our model all neurons are identical, and none can oscil- 
late autonomously, a few neurons, selected through the 
network connectivity, are timing the oscillations. This 
subpopulation of spike leaders can be interpreted as the 
emergent pacemakers of the network. The other neurons 
effectively amplify their action. 

Network degradation, i.e., randomly knocking out neu- 
rons, leads to complex changes in the set of these emer- 
gent pacemaker neurons. It also takes longer for them to 
reach threshold as N is reduced so the oscillation period 
increases. For lower values of AV{C — 0), i.e., weakly 
excitable networks, the period diverges along the phase- 
boundary between the SO and Q phases in Fig.[T}3, which 
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FIG. 3: (color online) Predicted rank: neurons ranked accord- 
ing to the number of neurons linked to them by one input link 
(green circles, "Level 1"), one or two input links (blue squares, 
"Level 2"), one, two or three input links (red diamonds, "Level 
3"). Actual rank: neurons ranked according to their somatic 
potential preceding a firing burst. Low rank corresponds to 
high potential. The neuron with lowest predicted rank indeed 
is the neuron with the lowest actual rank, but predicted rank 
correlates best with actual rank for higher rankings 




FIG. 4: (color online) fc-cores of a symmetric TV x TV random 
adjacency matrix. Nodes making up the 5-core are marked 
in red, 4-core nodes in blue, 3-core nodes in green, and 2-core 
nodes in orange. The radial distance of a node from the center 
increases with decreasing k. With a given fc-core, a node's 
radial position is increased by connections to lower fc-core 
nodes and decreased by connections to higher fc-core nodes. 
The four figures show a progressively increasing network size: 
(TVa = 40, TVs = 41, TVc = 42, No = 43). Image created 
using Ref . [13] 



agrees with the predictions of mean-field theory. For 
higher values of AV{C = 0) the system enters the HA 
phase. Unlike in mean-field theory, the HA phase of the 
random network exhibits the complex dynamical behav- 
ior, with period doubling and deterministic chaos, that 
was reported experimentally [B] . One example of this is 
shown in Fig. [2|3 where groups of high-activity bursts 
alternate with periods of deterministic chaos, forming a 
complex limit cycle. 

The SO <-> HA threshold curve TVc(AV") has a surpris- 
ing stair-case dependence on AV differing dramatically 
from the continuous curve predicted by the mean-field 
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theory (see Fig. [T]). These discontinuities define certain 
privileged numbers of neurons iVfe for which the network 
fails to support stable oscillations. The values of these 
iVfc's are independent of system parameters such as AC. 
The boundary of the SO regime makes discrete jumps be- 
tween Nk's as the parameters (AC, AV^) of the neuronal 
dynamical model are changed. In contrast, the phase 
boundary separating the SO and Q phases in Fig. [T]3 
largely follows mean-field predictions. 

Both the selection of the pacemaker neurons and the 
values of the privileged number of neurons appear to be 
determined largely by the mathematical properties of the 
adjacency matrix Mij independent of the details of our 
dynamical model. In Fig. |3] we show how one can iden- 
tify the pacemaker neurons through their connectivity by 
ranking them according to the number of neurons con- 
nected to a given neuron by no more than three links. 
These pacemaker neurons do not, however, play a cen- 
tral role in the determination of the minimum number 
of network neurons able to support collective oscillations 
for highly excitable neuronal networks. While sufficiently 
active synaptic inputs to a pacemaker indeed will more 
quickly drive the dendritic calcium concentration past 
the threshold C* , so that it becomes desensitized and 
thus unable to spike, a small set of such desensitized and 
quiescent neurons will not drive an active network to col- 
lectively desensitize. Rather, one must have a system- 
spanning high-connectivity network capable of simulta- 
neously desensitizing all of the neurons to quiet the in- 
herently excitable system. 

To quantify the size of a high-connectivity part of the 
network, it is useful to introduce the concept of a k- 
core [16]. A fc-core of a graph (for integer /c) is a sub- 
graph in which all nodes (i.e. neurons) have at least k 
inputs from other nodes in the subgraph. As the number 
of nodes increases in an Erdos-Renyi random network, 
/c-core clusters appear with larger k values at sharply 
defined thresholds. As an example we show in Fig. |4|\_ 
the fc-cores of a symmetric 40 x 40 random adjacency 
matrix. Nearly all nodes form a single 4-core cluster. 
Adding one more node at random does not change this 
feature (Fig. |4^), but adding two nodes at random, so 
that N = = 42 produces a sharp transition in which 
the network is now dominated by a single, system-sized 
5-core cluster (Fig. |4p). Adding an additional node to 
A'' = 43 does not alter the dominance of the 5-core, as 



shown by Fig. [4|D. For the random network used to gen- 
erate Fig. [l}3, these discontinuous transitions take place 
at A'3 = 17, when a 3-core appears; at iV4 = 26 when a 
4-core appears; at iVs = 37 when a 5-corc appears, and 
so on. The values of Nk for this realization of the random 
network are represented by dotted lines in Fig. [lj3. The 
locations of the discontinuities of the phase boundary as 
a function of N agree well, though not perfectly, with 
the fc-core transition values Nk- The discrepancies are 
presumably due to the fact that a member of a fc-core 
can have more than k input links, including links from 
non fc-core neurons. We emphasize that the fc-core con- 
cept is inapplicable to the SO <-*■ Q transition. Along the 
transition line, the few neurons with the highest connec- 
tivity are able to trigger an excitation wave that spreads 
through the whole system. These few emergent pacemak- 
ers are simply outliers having maximal connectivity and 
need not be part of a high fc fc-core. 

In summary, we have presented a simple model for 
rhythmogenic neuronal networks, such as the pBC, us- 
ing a combination of excitable integrate-and-fire neurons 
modified by a slower process of calcium-mediated desen- 
sitization. The most important conclusion of our work 
is that key features of the network dynamics - determi- 
nation of the pacemaker neurons and determination of 
the minimal number of neurons that support stable os- 
cillation - are determined (largely) by network connec- 
tivity. We also showed that in the phase diagram there 
is an asymmetry between the transition from the stable 
oscillation phase to the quiescent phase and the transi- 
tion from the stable oscillation phase to the high activity 
phase. The first transition is well described by mean- 
field theory, while the staircase structure of the phase 
boundary of the second transition refiects the full nature 
of network connectivity. This asymmetry originates from 
the difference between the dynamics of a spreading wave 
of voltage-mediated excitation and collective, calcium- 
mediated desensitization. Tests of the model should be 
straightforward. The excitability of neurons can be in- 
creased in experiment, effectively controlling the size of 
the action potential AV in our model. Measuring the 
onset action potential for complex firing patterns as a 
function of the number N of neurons should then directly 
reveal the predicted staircase structure of Fig. [l]3. 

We thank J. Feldman for enjoyable conservations and 
for sharing unpublished data on pBC dynamics. 
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